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Abstract 

In this paper we introduce a novel particle filter scheme for a class of partially- 
observed multivariate diffusions. We consider a variety of observation schemes, in- 
cluding diffusion observed with error, observation of a subset of the components of 
the multivariate diffusion and arrival times of a Poisson process whose intensity is a 
known function of the diffusion (Cox process). Unlike currently available methods, 
our particle filters do not require approximations of the transition and /or the obser- 
vation density using time-discretisations. Instead, they build on recent methodology 
for the exact simulation of the dif fusion process and th e unbiased estimation of the 
transition density as described in iBeskos et all l|2006tJ ). We introduce the Gener 



alised P oisson Estimator, which generalises the Poisson Estimator of IBeskos et al 
(|2006bh . A central limit theorem is given for our particle filter scheme. 



Keywords : Continuous-time particle filtering, Exact Algorithm, Auxiliary Variables, Cen- 
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1 Introduction 

There is considerable interest in using diffusion processes to model continuous-time phe- 
nomena in many diverse scientific disciplines. These processes can be used to model directly 
the observed data and/or to describe unobserved processes in a hierarchical model. This 
paper focuses on estimating the path of the diffusion given partial information about it. 
We develop novel particle niters for analysing a class of multivariate diffusions which are 
partially observed at a set of discrete time-points. 
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Particle filtering methods are standar d Monte-Carl o meth ods for analysing partially- 
observed discrete-time dynamic models (IDoucet et all 120011 ). They involve estimating 



the filtering densities of interest by a swarm of weighted particles. The approximation 
error decreases as the number of particles, N, increases. However, filtering for diffusion 
processes is significantly harder than for discrete-time Markov models since the transition 
density of the diffusion is unavailable in all but a few special cases. In many contexts even 
the observation density is intractable. Therefore, the standard propagation/weighting/re- 
sampling steps in the particle filter algorithm cannot be routinely applied. 

To circumvent these complications, a further appro ximation, based on a time-discretisat ion 



of th e diffusion, has been suggested (see for example ICrisan et all Il999l : iDel Moral et al. 



200ll ). The propagation of each particle from one observation time to the next is done by 



splitting the time increment into M, say, pieces and performing M intermediate simula- 
tions according to an appropriate Gaussian distribution. As M gets large this Gaussian 
approximation converges to the true diffusion dynamics. In this framework the compu- 
tational cost of the algorithm is of order M x N, and the true filtering distributions are 
obtained as both M and iV increase. 

Our approach does not rely on time-discretis ation, but builds on recent work on the 



Exact Algorit hm for the simulation of diffusions (IBeskos and Roberts! . 120051 ; iBeskos et al 



20063. l2005bl) and on the unbiased estimation of the diffusion transition density (IBeskos et al 



2006bl . l2005al ). This algorithm can be used in a variety of ways to avoid time discretisa- 



tions in the filtering problem. The potential of th e Exact Algorithm in the filtering problem 
was brought up in the discussion of IBeskos et al.l (l2006bl ). see the contributions by Chopin, 
Kiinsch, and in particular Rousset and Doucet who also suggest the use of a random weight 
particle filter in this context. 

One possibility is s imply to use the Ex act Algorithm to propagate the particles in the 
implementation of the I Gordon et al.l (119931 ) bootstrap particle filter, thus avoiding entirely 
the M intermediate approximate simulations between each pair of observation times. We 
call this the Exact Propagation Particle Filter (EPPF). Where possible, a better approach 
is to adapt the Exact Algorithm to simulate directly from (a particle approximation to) 
the filtering density using rejection sampling; we term this the Exact Simulation Particle 
Filter (ESPF). 

However, our favoured meth od goes in a different direc tion. We work in the framework 
of the auxiliary particle filter of iPitt and Shephardl (ll999l ). where particles are propagated 
from each observation time to the next according to a user-specified density and then are 
appropriately weighted to provide a consistent estimator of the new filtering distribution. 
Due to the transition density being unavailable, the weights associated with each particle 
are intractable. However, our approach is to assign to each particle a random positive 
weight which is an unbiased estimator of the true weight. We call this the Random Weight 
Particle Filter (RWPF). Our algorithm yields consistent estimates of the filtering distribu- 
tions. The replacement of the weights in a particle filter by positive unbiased estimators is 
an interesting possibility in more general contexts than the one considered in this paper. 
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Indeed, in Section 13.21 we show that this approach amounts to a convenient augmentation 
of the state with auxiliary variables. 

The construction of the unbiased estimators of the weights is one of the main contri- 
butions of this paper, and it is of indep e ndent interest. This is based on an extension 
of the Poisson Estimator of iBeskos et al.l (j2006bl ). which we call the Generalised Poisson 
Estimator. This estimator is guaranteed to return positive estimates (unlike the Poisson 
Estimator) and its efficiency (in terms of variance and computational cost) can be up to 
orders of magnitude better than the Poisson Estimator. Optimal implementation of the 
Poisson and the Generalised Poisson estimators is thoroughly investigated theoretically 
and via simulation. 

All three time-discretisation- free particle filters we introduce are easy to implement, 
with the RWPF being the easiest and the most flexible to adapt to contexts more general 
than those considered here. A simulation study is carried out which shows that the RWPF 
is considerably more efficient than the ESPF which is more efficient than the EPPF. We also 
provide a theoretical result which shows that our filters can have significant computational 
advantages over time-discretisation methods. We establish a Central Limit Theorem (CLT) 
for the estimation of expectations of the filtering distributions using e it her o f the EPPF, 
ESPF and the RWPF. This is an extension of the results of IChopinl (120041 ) . The CLT 
shows that, for a fixed computational cost K, the errors in the particle approximation of 
the filtering distributions decrease as K~ x l 2 in our methods, whereas it is known that the 
rate is K~ 1 ^ or slower in time-discretisation methods. 

The main limitation of the methodology presented here is the requirement that the 
stochastic differential equation specifying the underlying diffusion process can be trans- 
formed to one with orthogonal diffusion matrix, and gradient drift. Although this frame- 
work excludes some important model types (such as stochastic volatility models) it incor- 
porates a wide range of processes which can model successfully many physical processes. 
On the other hand, our methods can handle a variety of discrete-time observation schemes. 
In this paper we consider three schemes: noisy observations of a diffusion process, observa- 
tion of a subset of the components of a multivariate diffusion, and arrival times of a Poisson 
process whose intensity is stochastic and it is given by a known function of a diffusion. 

The paper is organised as follows. Section [2] introduces the model for the underlying 
diffusion and the necessary notation, the observation schemes we consider and the simulated 
data sets on which we test our proposed methods. Section [3] introduces the RWPF and 
states the CLT. Section H] introduces the main tool required in constructing the RWPF, 
the Generalised Poisson Estimator (GPE). Several theoretical results are established for 
the GPE, and a simulation study is performed to assess its performance. Section [5] is 
devoted in the empirical investigation of the performance of the different particle filters 
we introduce. Several implementation issues are also discussed. Section [U] closes with a 
discussion on extensions of the methodology and the appendices contain technical results 
and proofs. 
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2 Signal, data and assumptions 

Let the signal be modelled by a (i-dimensional diffusion process 

dX s = a(X fl ) ds + dB s , s G [0, t] . 



(1) 



We assume throughout the paper that the drift is known. Our approach requires some 
assumptions which we summarize in this paragraph: i) a is continuously differentiable in 
all its arguments, ii) there exists a function A : R d — > R such that a(u) = VA(u), and 
hi) there exists / > — oo such that 0(u) := (||q:(u)|| 2 + V 2 A(u))/2 — I > 0. Among these 
last three conditions i) and iii) are weak and the strictest is ii), which in the ergodic case 
corresponds to X being a time-reversible diffusion. 

The transitio n density of (ED) is typically intractable but a useful expression is av ailable 



(see for example iBeskos et al. 



[VP 

2006bl ; iDacunha-Castelle and Florens-Zmiroul . Il986f ) 



Pt(x t | x ) = A/" t (x< - x ) exp{A(x t ) - A(x ) - lt}E 



exp 



0(W s )ds 



(2) 



In this expression A/t(u) denotes the density of the d- dimensional normal distribution with 
mean and variance tlj evaluated at u G R d , and the expectation is taken w.r.t. a 
Brownian bridge, W s , s G [0, t], with Wo = xo and W t = x t . Note that the expectation 
in this formula typically cannot be evaluated. 

The data consist of partial observations y\, y 2 , . . . , y n , at discrete time-points < t\ < 
■ ■ <t n . We consider three possible observation regimes: 

Diffusion observed with error. The observation y^ is related to the signal at time £j 
via a known density function /(y^x^). This model extends the general state-space 
model by allowing the signal to evolve continu ously in time. Ther e is a wide range 
of applications which fit in this framework, see iDoucet et al.l (120011 ) for references. 



to < 



(A) 



(B) 



(C) 



Partial Information. At time ti we observe yi = CP^-u) f° r some non- invert ible known 
function ((■). For example we may observe a single component of the d- dimensional 
diffusion. In this model type /(yi|xtj = 1 for all x ti for which C( X *J — V%- 

Cox Process. In this regime the data consist of the observation times ti which are 
random and are assumed to be the arrivals of a Poisson process of rate z/fX,,), for 
some known functi on v. Such models are popular in insur ance (IDassios and Jang . 



20051 1 and finance (jEngell . l2000l : iDuffie and Singleton! . Il999l ). and they have r ecently 



been used to analyse data from single molecule experiments (IKou et all 120051 ) .There 
is a significant difference between this observation regime and the two previous ones. 
To have notation consistent with (A) and (B) we let y^ = t t denote the time of the ith 
observation; and define the likelihood f(yi \ x^ i _ 1 ,x ti ) to be the probability density 
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that the next observation after ti-i is at time £j. This density can be obtained by 
integrating 

KXJexpj-jT KX s )d S |, (3) 

w.r.t. the distribution of (X a ,s G conditionally on X fi _ 1 = x^ ^Xt. = x fi . 

The distribution of this conditioned process has a known dens i ty w.r. t. the Brownian 



bridge measure and it is given in Lemma 1 of iBeskos et al.l (l2006bl ). We can thus 
show that the density of interest is 



, : - exp{A(x t J-A(x ti _ 1 )}E 

Pu-ti-A^ti I x tj _J 



expj-^ (0(W S ) + !/(W s ))ds 



'ti-l 

(4) 

where expectation is with respect to the law of a Brownian Bridge from x ti _j to x^. 

We take a Bayesian approach, and assume a prior distribution for X . Our interest 
lies in the online calculation of the filtering densities, the posterior densities of the signal 
at time ti given the observations up to time £j, for each 1 < i < n. While these densities 
are intractable, we propose a particle filter scheme to estimate recursively these densities 
at each observation time-point. As we point out in Section 6, our approach allows the 
estimation of the filtering distribution of the continuous time path (X s ,tj_! < s < tj). 

A more flexible model for the signal is a diffusion process Z which solves a more general 
SDE than the one we have assumed in (CD): 

dZ s = b(Z s )ds + S(Z,)dB s , se[0,t]. (5) 

In contrast with (CD), ([5]) allows the diffusion coefficient to be state-dependent. Our methods 
directly apply to all such processes provided there is an explicit transformation Z s i— ► 
tj(Z s ) =: X s , where X solves an SDE of the type ([T]); the implied drift ct can be easily 
expressed in terms of b and I] via Ito's formula and it will have to satisfy the conditions we 
have already specified. In model (A) the likelihood becomes f(yi \ r\ 1 (X ti )) , in model (B) 
the data are yi = C(?7^ 1 (X t J) and in model (C) the Poisson intensity is z/(?7 _1 (X s )), where 
rj 1 denotes the inverse transformation. Therefore, the extension of our methodology 
to general diffusions is straightforward when d — 1; under mild conditions © can be 
transformed to ([T]) by r](Z s ) = J^ B d^, for some arbitrary u* in the state space of 

the diffusion. Moreover, the drift of the transformed process will typically satisfy the three 
conditions we have specified. However, the extension is harder in higher dimensions. The 
necessary transf ormation is more c omplicated when d > 1 and it might be intractable or 



even impossible (jA'it-Sahalial . 12004 ). Even when such a transformation is explicit it might 
imply a drift for X which violates condition ii). Nevertheless, many physical systems can 
be successfully modeled with diffusions which can be transformed to (DO). 

Our particle filtering methods will be illustrated on two sets of simulated data: 
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Example 1: Sine diffusion observed with error. The signal satisfies 



dX s = sm(X s )ds + dB s , (6) 

and the data consist of noisy observations, ~ N(X^,cr 2 ). Figure [TJ^top) shows a simula- 
tion of this model with a = 0.2. In this case 

(f)(u) = (sin(w) 2 + cos(w) + l)/2 . (7) 

This process is closely related to Brownian motion on a circle. It is convenient as an 
illustrative example since discrete-time skeletons can be ea sily simulated from th is process 



using the most basic form of the Exact Algorithm (EAl in lBeskos et all l2006al . R-code is 
available on request by the authors). 

Example 2: OU-driven Cox Process. The second data set consists of the arrival times 
of a Poisson process, y^ = ti, whose intensity is given by v{X s ), s > 0, where 

v{x) = a + P\x\, 

and X is an Ornstein-Uhlenbeck (OU) process, 

dX s = -pX s ds + dB s . 

The OU process is stationary with Gaussian marginal distribution, N(0, 1/(2/?)). Thus, 
an interpretation for this model is that the excursions of X increase the Poisson intensity, 
whereas a corresponds to the intensity when X is at its mean level. An example data set is 
shown in Figure EJ where we have taken a = 0, (3 = 20, p = 1/2. Although the transition 
density of the OU process is well-known, 



X t | X = x ~ N ( e'Ptxo, ^(1 - e^] 



the observation density f(yi+i \ x ti ,x ti+1 ) is intractable. 

Examples 1 and 2 are examples of observation regimes (A) and (C) respectively. We 
will show that observation regime (B) can be handled in a similar fashion as (A), so we 
have not included an accompanying example. 



3 Random weight particle filter 

As in Section [2] we will denote the observation at time tj by y^, and p t (- \ ■) will denote the 
system transition density over time t (see Equation^. We will write A$ = t i+ i — t iy and 
the filtering densities p(x t Jy 1: j) will be denoted by 7Tj(x t .), where by standard convention 
Vv.i = (yi, ■ ■ ■ ,Vi)- To simplify notation, when we introduce weighted particles below, we 
will subscript both particles and weights by i rather than ti. 
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Our aim is to recursively calculate the filtering densities 7i"i(x ti ). Basic probability 
calculations yield the following standard filtering recursion for these densities 

7r m (x t . +1 )oc J /(i/i+xIxt^x^JpA^x^JxtjTr^XtJdx^. (8) 

Particle filters approximate 7Ti(x ti ) by a discrete distribution, denoted by 7Tj(xtJ, whose sup- 
port is a set of iV particles, {x^}^, with associated probability weight {w^}jL v Substi- 
tuting VTj(x t .) for 7Tj(x t .) in (JH]), yields a (continuous density) approximation to 7r i+1 (x ti+1 ), 



N 



7T 



Si) 



(9) 



The aim of one iteration of the particle filter algorithm is to construct a further particle 
(discrete distribution) approximation to 7fj + i (x ti+1 ) . 

We can obtain such a particle approximation via importance samp ling, and a general 
frame work for achieving this is given by the auxiliary particle filter of iPitt and Shephard 
( 119991 ). We choose a proposal density of the form 



N 



(10) 



Choice of suitable proposals, i.e. choice of the fl^'s and q, is discussed in the analysis of 
our specific applications in Section [51 

To simulate a new particle at time t i+ i we (a) simulate a particle x^ at time i, where 
k is a realisation of a discrete random variable which takes the value j G {1,2,..., N} with 



probability , and (b) simulate a new particle at time t i+ % from g(x ti+1 |x^ j , y i+1 ). The 
weight assigned to this pair of particles (x^ , x fi+1 ) is proportional to 



w? f (Vi+i I , x ti+1 )p Ai (x t i+1 1 x 



^ k) q{x ti+1 \x.\ K> ,y i+1 



(fc) 



(11) 



rO') ,.,0') ^ 



This is repeated N times to produce the set of weighted particles at time 

which gives an importance sampling approximation to 7r i+1 (x^ +1 ). Renormalising the 
weights is possible but does not materially affect the methodology or its accuracy. Im- 
provem ents on independent sa mpling in step (a) can be made: see the stratified sampling 
ideas of lCarpenter et al.l (119991 ) . The resulting particle filter has good theoretical properties 
including consistency (ICrisanl. 120011) and central limit theorem s for e stimates of posterior 
moments foel Moral and Miclol . l2000l : IChopinl . I2004J : iKiinschl . l2005h . as N -> 00. Un- 
der conditions relati ng to exponential forgett i ng of initial conditio ns, particle filter errors 
stabilise as n — > 00 ( iDel Moral and Guionnetl . l200ll ; IKiinschl . 120051 ). 



N 
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The difficulty with implementing such a particle filter when the signal X is a diffusion 
process is that the transition density PAi( x t i+1 ) which appears in ffTTl) is intractable for 
most diffusions of interest, due to the expectation term in ([2]). Furthermore, for observation 
model (C) (but also for more general models), the likelihood term f(yi + i\x\ k \ x ti+1 ) given 
in (TJ1) cannot be calculated analytically. 

We circumvent these problems by assigning each new particle a random weight which is 
a realisation of a random variable whose mean is (fTTl) . The construction and simulation of 
this random variable is developed in Section HI and it is based on the particular expression 
for the transition density in §2§. The replacement of the weights by positive unbiased 
estimators is an interesting possibility in more general contexts than the one considered 
in this paper. Indeed, in Section 13.21 we show that this approach amounts to a convenient 
augmentation of the state with auxiliary variables. 

3.1 Simulation of weights 

In all models the weight associated with the pair (x- fc \x 4i+1 ) equals 

^i+i(xf \ x* j+1 , y i+ i)/i fl (xf } , x t . +1 , U, t i+ i) (12) 
where h i+ i is a known function, and for < u < t, 



H g (x,z,u,t) := E 



exp{-jf <?(W s )d S J 



where the expectation is taken w.r.t. a <i-dimensional Brownian bridge W, starting at time 
u from W M = x and finishing at time t at W t = z. 



Models (A) and (B) : For these model types 

, / (fc) \ Wi k) f{yi+i\xt 1+1 )AfAA^ . 



^V(?/m|xt I+1 )AUx fl+1 - xg) exp{A(x, +1 ) - A(xf } )} 
Pi <H x t i+ il x i iVi+i) 



i + l \ I 

and g = (p. In model type (B) the proposal distribution g(x ti+1 |x- fc \ yi + i) should be chosen 
to propose only values of x ti+1 such that C( x t l+1 ) = Vi+i] then f(y i+ x\x ti+1 ) = 1. 

Model (C): A synthesis of ©, Q and ([II]), with g = <\> + v gives 

, / (*) x _ wf )z/ (x* I+1 )A/'A I (x tl+1 -xf ) )exp{A(xi i+1 ) - A(xf } )} 

h l+1 {^ ,x ti+1 ,y i+1 j - , 

Pi 9ix ti+1 |x i ,y l+ i) 

Section H] shows how to construct for each pair of (x, z) and times (u,t), with u < t, ad- 
ditional auxiliary variables V, and a function r(v,x,z, u, t) > 0, with the property that 
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E[r(V, x, z, u, t) | x, z] = /i s (x, z, u, t). The auxiliary variables are simulated according to 
an appropriate conditional distribution Q g (- | x, z, u, t), and r is easy to evaluate. Our 
method replaces in the weight the intractable term fi g with its unbiased estimator r. 

Random Weight Particle Filter (RWPF) 

PFO Simulate a sample . . . , xj, from p(x ), and set = 1/N. 
For i — 0, . . . , n — 1, for j = 1, . . . , iV: 

PF1 calculate the effective sample size of the {$* ] }, ESS = (Ej^Oflf ) 2 ) -1 ; i/^SS < C, 
for some fixed constant C, simulate k iy j from p(k) = (3\ k \ k = 1, ...,N and set 
8^ = 1; otherwise set kij = j and 5^ = pf ; 

PF2 simulate x^_\ from ?(x ti+1 yi+i); 

PF3 simulate v i+1 ~ Q g ( - \ xf x t . +1) t i? 

PF^ assign particle x^ a weight 

w&i = 4+ } i^+i(xf xg 1; y i+ i)r(v i+1 , xf w) , x ti+1 , t i} . (13) 



Notice that this algorithm contains a decision as to whether or not resample particles 
before propagation in step PF1, with decision being based on the ESS of the (3\ j) . The 
constant C can be interpreted as the minimum acceptable effective sample size. (See Liu 
and Chen 1998 for the rationale of basing resampling on such a condition.) Whether or 
not resampling occurs will affect the weight given to the new sets of particles, and this is 
accounted for by different values of 8jV\ in PF1. Optimally, the resampling for step PF1 
will incor porate dependen c e acr oss the N samples; for exa mple the strat i fied s ampling 
scheme of ICarpenter et al.l (119991 ) or the residual sampling of iLiu and Chenl (119981 ). 



3.2 An equivalent formulation via an augmentation of the state 

In the previous section we described a generic sequential Monte Carlo scheme where the 
exact weights in the importance sampling approximation of the filtering distributions are 
replaced by positive unbiased estimators. We now show that this scheme is equivalent to 
applying an ordinary auxiliary particle filter to a model with richer latent structure. We 
demonstrate this equivalent representation for model types (A) and (B), since an obvious 
modification of the argument establishes the equivalence for model type (C). 

According to our construction, conditionally on X^, X ti+1 , U and t i+ i, VVfj is inde- 
pendent of Vj and Xj . for any j different from + Additionally, it follows easily from 
the unbiasedness and positivity of r that, conditionally on X tj = x, r(v i+1 , x, x ti+1 , t i: t i+ i) 
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is a probability density function for (X ti , Vj+i) with respect to the product measure 
Lefe(dz) x Q g (dv | x, z, tj, where Left denotes the Lebesgue measure. 

Consider now an alternative discrete-time model with unobserved states (Zj, Vj),i = 
1, . . . , n, Zj G R d , with a non-homogeneous Markov transition density 

p i+1 (z i+1 , V i+1 | Zj,Vj) = r(v i+1 ,Zi,Zi +1 ,ti, t i+ i) , 

(this density is with respect to Leb x Q g ) and observed data yi with observation density 
f{yi+i | Zj,Zj + i). By construction the marginal filtering distributions of Zj in this model 
are precisely 7Tj(x t .), i.e. the filtering densities in Consider an auxiliary particle filter 
applied to this model where we choose with probability each of the existing particles 
(zp , v^), and generate new particles according to the following proposal 

(z i+ i,v i+ i) ~ q(z i+1 | z ( i k \y i+1 )Q g (d\ i+ i \ zf \ z i+1 , U, t i+l )Leb(dz i+1 ) , 

where q is the same proposal density as in f TTOT) . The weights associated with each particle 
in this discrete-time model are tractable and are given by ( |T3l) . Therefore, the weighted 

sample < (z^_\, ^i+i) [ is precisely a particle approximation to 7Ti + i(x ti+1 ), and RWPF is 

equivalent to an auxiliary particle filter on this discrete-time model whose latent structure 
has been augmented with the auxiliary variables V*. 

This equivalent representation sheds light on many aspects of our method. Firstly, it 
makes it obvious that it is inefficient to average more than one realization of the positive 
unbiased estimator of \i g per particle. Instead it is more efficient to generate more particles 
with only one realization of the estimator simulated for each pair of particles. 

Secondly, it illustrates that RWPF combines the advantages of the bootstrap and the 
auxiliary particle filter. Although it is easy to simulate from the probability distribution 
Q g (as described in Section HJ), it is very difficult to derive its density, (with respect to an 
appropriate reference measure). Since the VjS are propagated according to this measure, its 
calculation is avoided. This is an appealing feature of the bootstrap filter which propagates 
particles without requiring analytically the system transition density. On the other hand 
the propagation of the ZjS is done via a user-specified density which incorporates the 
information in the data. 

Thirdly, it suggests that the RWPF will have similar theoretical properties with auxil- 
iary particle filters applied to discrete-time models. This is explored in Section 13.31 

3.3 Theoretical properties 

Consider estimation of the posterior mean of some function (p of the state at time ti, 
E[y(xt i )|j/i : j]. A natural approach to the investigation of particle filter effectiveness is to 
consider the limiting behaviour of the algorithm as N — > oo. For the standard auxiliary 
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particle filter, IChopinl (120041 ) introduces a central limit theorem (CLT) for estimation of 
this type of expectations. This CLT applies directly t o bot h EPPF and the ESPF. 

In Appendix F we extend the result of IChopinl (120041 1 and give a further necessary 
condition on the random weights in RWPF under which a CLT still holds. This extra 
condition is (C2). The expression for the variance of the estimator of E[<^(xt i )|yi : j] obtained 
with RWPF differs from the expression in the standard case (i.e. when the weights are 
known) by an extra term caused by the randomness in the weights (see Equations I2TH2TH 
and the comment on Theorem [3] in Appendix F for further details) . The ready adaptation 
of Chopin's approach is facilitated by the observation that the RWPF can be re-expressed 
as a standard particle filter for the an augmented state (see Section 13.21) . 

One important consequence of this CLT is that the errors in estimating E[yj(x t J|?/ 1: j] 
are of order iV -1 / 2 . Previous filtering methods for the diffusion problems we consider are 
based on (i) discretising time and introducing M intermediate time points between each 
observation time; (ii) using an Euler, or higher order, approximation to the diffusion (CQ); 
and (iii) ap plying a particle, or other, filter to this approximate discrete time model. Sec 
for example 
are given by 



Crisan et al.l (119991) . R esults giving the order of the errors in one such scheme 



Del Moral et al.l (120011 ). For Models such as (A) and (B) the errors are of order 



N x l 2 provided that the number of intermediate time steps M between each observation 
increases at a rate N 1 / 2 . Thus for fixed computational cost K oc MN the errors decrease 
at a rate K~ l l^. For models such as (C), where the likelihood depends on the path of the 
state between two s uccessive observations, the rate at which errors de crease will be slower , 
for example K' 1 ^ dPel Moral et all l200lh . or K' 1 / 6 (Theorem 1.1 of lCrisan etaflll999h . 



4 Generalised Poisson Estimators 

We have already motivated the need for the simulation of a positive unbiased estimator of 

E [E] where E =: exp j- jf #(W s )ds j , (14) 

where the expectation is taken w.r.t. a <i-dimensional Brownian bridge W. In this section 
we introduce a methodology for deriving such estimators, and provide theoretical and 
simulation results regarding the variance of the suggested estimators. These results are of 
independent interest beyond particle filtering, so we present our methodology in a general 
way, where g is an arbitrary function assumed only to be continuous on H d . We assume 
that Wo = x and W$ = z, for arbitrary x, z G R d and t > 0. By the time- homogeneity 
property of the Brownian bridge our methodology extends to the case where the integration 
lim its change to u and u + 1, for any u > 0. 



Beskos et al.l (j2006bl ) proposed an unbiased estimator of (1T41) . the Poisson Estimator. 



K 

(Mt A- K n[ c -5( w j] ; ( 15 ) 



PE 

3=1 
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K is a Poisson random variable with mean At, the ipjS are uniformly distributed on [0, t], 
and c G R, A > are arbitrary constants. (Here and below we assume that the empty 
product, i.e. when k = 0, takes the value 1.) The two main weaknesses of the PE are 
that it may return negative estimates and that its variance is not guaranteed to be finite. 
Both of these problems are alleviated when g is bounded. However this is a very restrictive 
assumption in our context. Therefore, here, we introduce a collection of unbiased and 
positive estimators of (1141) which generalise the PE. The methods we consider allow c and 
A depend on W, and permit k to have a general discrete distribution. Firstly, we need to 
be able to simulate random variables Lw and U-w with 



£w < #(W S ) < U w , for all s G [0, t], 



(16) 



and to be able to simulate W s at any s, given the condition implied by (U6l) . For un- 
bounded g this is non-trivial. However, both of these simula t ions ha ve become feasible 
since the introduction of an efficient algorithm in iBeskos et al. fl2005bh . An outline of the 
construction is given in Appendix A. 

Let U-w and L w satisfy ffTB"]) and ipj,j > 1, be a sequence of independent uniform 
random variables on [0, t]. Then, ( TT4"1) can be re-expressed as follows, 



E 



e ^^exp 



(U w -g{W s ))ds 



E 



E 



k=0 



t \ k ' 

(U w -g(W s ))ds 



k=0 ' j=l 



E 



l[(U w -g(W^)) 



3=1 



;i?) 



where k is a discrete random variable with conditional probabilities P[k = k \ Uw, L-w] — 
p(k | Uw,Lw). The second equality in the above argument is obtained using dominated 
convergence and Fubini's theorem (which hold by positivity of the summands). 

We can derive various estimators of (1T4")) by specifying p(- \ U-w, L w ). The family of all 
such estimators will be called the Generalised Poisson Estimator (GPE): 



GPE: e- Uwt 



t K 



k\p(k I U w , L w 



l[(U w - g(W^)) . 



The following Theorem (proved in Appendix B) gives the optimal choice for p(- | U-w, L-w) 
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Theorem 1. The conditional second moment of the Generalised Poisson Estimator given 
£/w an d £w; is: 



-2Uwt 



E 

k=0 



p(k | f/w, L^)k\ 2 



E 



j\u w -g(W s )) 2 ds^ \U Wl L w 



(19) 



If 



« t k/2 



fc=0 



t N * 

2. 



(f/w-^W,))^ |C/ W ,L W 



1/2 



< oo 



(20) 



i/ien t/ie second moment is minimised by the choice 



t k/2 

p(k | U W ,L W ) oc — E 



j\u w -g(W s )) 2 ds^J \U W ,L W 



1/2 



(21) 



with minimum second moment given by 

t k/2 



-Uwt 



k=0 



r* \ k 1 1/2N 

j o (U w -g(W s )) 2 dsj \U Wl L w 



< oo , /or almost all U-w, L w . 

(22) 



Whilst the right-hand side of ( J2TT) cannot be evaluated analytically, it can guide a suitable 
choice of p(- \ Uw, Lw). If W were known, the optimal proposal is Poisson with mean 



A w := (t JjU w -g(W s )) 2 ds^ 



1/2 



(23) 



We will discuss two possible ways that (1231) can be used to choose a good proposal. 

A conservative approach takes p(- \ Uw, Lw) to be Poisson with mean (C/w — -^w)^ (an 
upper bound of Aw)- We call this estimator GPE-1. An advantage of GPE-1 is that its 
second moment is bounded above by E[e _2iw *]. Thus, under mild and explicit conditions 
on g, which are contained in the following theorem (proved in Appendix C), the variance 
of the estimator is guaranteed to be finite. 

Theorem 2. A sufficient condition for GPE-1 to have finite variance is that 



g(ui, ...,u d/ 



> -5 J^(l + \Ui\), for all Ui G R, 1 <i < d, 5>0. 



i=l 
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Since Aw is stochastic, an alternative approach is to introduce a (exogenous) random 
mean and assume that p(-\Uw, Lw) is Poisson with this random mean. For tractability we 
choose the random mean to have a Gamma distribution, when p(- | Uw, Lw) becomes a 
negative-binomial distribution: 

GpE . 2; e^S^nk - 0(wj|, (24) 

where 7 W and (3 denote the mean and the dispersion parameter respectively of the negative 
binomial. Since the negative-binomial has heavier tails than the Poisson Estimator, GPE- 
2 will have finite variance whenever there exists a PE with finite variance. On the other 
hand, big efficiency gains can be achieved if 7 W is chosen to be approximately E[Aw I 
U-WjL-w]. There is a variety of ad-hoc methods which can provide a rough estimation of 
this expectation. Applying Jensen's inequality to exchange the integration with the square 
power in (J23l) . and subsequently approximating E[g(W s ) | Uw, Lw] by p(E[W s ]), suggests 
taking 

7 W = tu w - g (*—f^ + yf J ds > ■ ( 25 ) 

A simulation study (part of which is presented in Section I4TT1 below) reveals that this choice 
works very well in practice and the GPE-2 has up to several orders of magnitude smaller 
variance than the PE or the GPE-1. The integral can usually be easily evaluated, otherwise 
a crude approximation can be used. 

We have confined our presentation to the case where the ex pectation in (fl4"|) is w .r.t. 
the Brownian bridge measure. Nevertheless, as pointed out in iBeskos et al. f|2006bl ) the 



PE can be constructed in exactly the same way when the expectation is taken w.r.t. an 
arbitrary diffusion bridge measure, as long as exact skeletons can be simulated from this 
measure. The GPE can also be implemented in this wider framework, provided that the 
process W can be constructed to satisfy (fTBjl . 

4.1 Simulation study 

We consider a smooth bounded test function g(u) = (sin(w) 2 -|-cos(u) + i)/2. This has been 
chosen in view of Example 1. The function g is periodic, with period 2tt. In [0,27r] it has 
local minima at and 2ir, global minimum at tc and maxima at 7r/3 and 5ir/3. Since g is 
bounded by 9/8 we can construct a PE which returns positive estimates by setting c > 9/8. 



Under this constraint, IBeskos et al.l (j2006bl ) argued that a good choice isc = A = 9/£ 



Simulation experiments suggested that the performance of the GPE-2 is quite robust to 
the choice of the dispersion parameter {3. We have fixed it in our examples to (3 = 10. Table 
[TJ summarizes estimates of the variance of the estimators based on 10 4 simulated values. 
We see that GPE-2 can be significantly more efficient than PE, in particular when taking 
into account E[/d. In general, the performance of PE is sensitive to the choice of c and A. 
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GPE-1 


4.21 x 10' 


-3 


0.208 


0.034 


GPE-2 


2.08 x 10~ 


-3 


0.220 


0.033 


Var(£) 


3.74 x 10" 


-5 


3.27 x 10~ 3 


4.72 x 10~ 3 


E[«] PE 


1.118 




1.126 


1.121 


GPE-1 


0.130 




1.091 


0.744 


GPE-2 


0.119 




0.329 


0.735 



Table 1: Monte Carlo estimates of the variance of four estimators of (fl4l) where giu) = 
(sin(w) 2 + cos(w) + l)/2. For comparison we give also var(.E'). We also report an estimate 
of E[ac]. We consider three different pairs of starting and ending points (x,z) and time 
increment t — 1. The estimates in the table were obtained from a sample of 10 4 realisations. 



GPE-1 is typically less efficient than GPE-2. Table 1 also gives the value of Vax(E) which 
takes significantly smaller values (by a couple of orders of magnitude) than any of PE, 
GPE-1 or GPE-2, illustrating the efficiency cost of these auxiliary variable constructions 
in absolute terms. 

We have also investigated how the efficiency of the PE and GPE-2 varies with the 
time increment t and in particular for small t (results not shown). These empirical results 
suggest that the coefficient of variation of the errors of both PE and GPE-2 are 0(t s ) for 
some 5 > 0; but that the value of 5 differs for the two estimators. In the cases that we 
investigated, the GPE-2 appears to have a faster rate of convergence than PE. 

The results of this simulation study have been verified for other functions g (results not 
shown). We have experimented with differentiable (e.g. g{u) = u) and non-differentiable 
(e.g. g(u) = \u\) unbounded functions. In these cases it is impossible to design a PE 
which returns positive estimates w.p.l. Again, we have found that the GPE-2 performs 
significantly better than the PE. 

It is important to mention that alternative Monte Carlo methods exist which yield 
consistent but biased estimates of fll4|) . One such estimator is obtained by replacing the 
time-integral in (JHj) with a Riemann approximation based on a number, M say, of interme- 



diate points. This technique is used to construct a transition density estimator in iNicolau 



and effectively underlies the transition density estimator of iDurham and Gallant 



2002) (when the diffusion process has constant diffusion coefficient). The approach of 



Durham and Ga 



20061 : 



Chib et al. 



lantl (|2002h has been used in MCMC and filtering applications ( IGolightly and Wilkinson! . 



20061 ; llonided . 120031 ) . In the filtering context it provides an alternative to 



RWPF, where the weights are approximated. It is not the purpose of this paper to carry 
out a careful comparison of RWPF with such variants. However, as an illustration we 
present a very small scale comparison in the context of estimating the transition density, 
Pt(z | a;), of (jSJ) for t — 1 and x, z as in Table [U We compare 4 methods. Two are based on 



15 



Estimator 


x = 0, z = 


X = 0, Z = 71 


X = TV, Z = 71 


PE 


1.25 


0.93 


0.17 


GPE-2 


0.13 


0.78 


0.2 


DG-1 


0.5 


0.45 


0.3 


DG-5 


0.28 


0.19 


0.22 



Table 2: Monte Carlo estimates based on 10 4 realisations of the root mean square error 
divided by the true value of 4 estimators of p t (z \ x), of (ED for t — 1 and various x, z. As 
true value we take the estimate produced by averaging the estimations given by GPE-2. 
The number of intermediate points used for each estimator are 1 and 5 for DG-1 and DG-5 
respectively; the number of Brownian bridge simulations for PE and GPE-2 are given in 
Table |H(E[k]). 



(j2J) and use the PE and the GPE-2 to generate es timators of the expectation. The other 
two, DG-1 and DG-5 are two implementation of the lDurham and Gallantl (120021 ) estimator, 
with 1 and 5 respectively intermediate points. We compare the methods in terms of their 
root mean square error divided by the true value (i.e. the coefficient of variation). As the 
true value we used the estimate of the GPE-2. The results of the comparison are presented 
in Table [2j Notice that DG-1 and DG-5 simulate many more variables than GPE-2 to 
construct their estimates. 



5 Comparison of particle filters on the simulated data 

We now demonstrate the performance of the different particle filters we have presented on 
the two examples introduced in Section [21 

5.1 Analysis of the sine diffusion 

We first consider analysing the sine diffusion of Example 1. The simulated data is shown 
in Figure [H(top). We compare four implementations of the particle filter each of which 
avoids time-discretisations by using methodology based on the Exact Algorithm (EA) for 
simulating diffusions: i) EPPF, which uses EA for implementing a bootstrap filter, ii) 
ESPF, which adapts EA to simulate by rejection sampling from the filtering densities, iii) 
RWPF1, an implementation of RWPF using PE (see Tabled]) to simulate the weights, iv) 
RWPF2, an implementation of RWPF using GPE-2 to simulate the weights. Details on 
the implementation of EPPF and ESPF are given in Appendix D. 

In this simple example ESPF is more efficient than EPPF, since it has the same compu- 
tational cost, but it is proposing from the optimal proposal distribution. However, we have 
efficiently implemented ESPF exploiting several niceties of this simple model, in particular 
the Gaussian likelihood and the fact that the drift is bounded. In more general models 
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Figure 1: Top: A realisation of the sine diffusion (black line) on [0, 100]; 100 observations 
at unit time intervals (blue circles); mean of the filtering distribution of the diffusion at the 
observation times obtained by RWPF2 with N = 1, 000 particles (red circles). Bottom: the 
difference between observed data and filtered means (circles) and 90% credible intervals 
(red dashed lines) from RWPF2. (Whilst for clarity they are shown for all times, the 
credible intervals were only calculated at the observation times.) 



implementation of ESPF can be considerably harder and its comparison with EPPF less 
favorable due to smaller acceptance probabilities. 

In this context where <j) is bounded one can speed up the implementation of GPE-2 with 
practically no loss of efficiency by replacing Uw in ([21]) and f[2"5"]) by 9/8 which is the upper 
bound of 0. In this case, there is no need to simulate t/w and Lw We have implemented 
this simplification in the RWPF2. 



Algorithms EPPF, RWPF1-2 used the stratified re-sampling algorithm of lCarpenter et al. 



( 119991 ). with re-sampling at every iteration. For RWPF1-2 we chose the proposal distribu- 
tion for the new particles based on the optimal proposal distribution obtained if the sine 
diffusion is approximated by the Ozaki discretisation scheme (details in Appendix E). For 

(k) 

EPPF we chose the f3\ s to be those obtained from this approximation. 

The number of particles used in each algorithm was set so that each filter had com- 
parable CPU cost, which resulted in 500, 500, 910 and 1000 particles used respectively 
for each algorithm. For these numbers of particles, EPPF and ESPF on average required 
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Figure 2: Relative efficiency of the 4 particle filter algorithms at estimating the filtering 
mean E[X ti \yi-.i]. Each line gives the relative efficiency of one algorithm compared to 
RWPF2 (black: RWPF2, green: RWPF1, red: ESPF, blue: EPPF). See text for details. 



the proposal of 1360 particles and required 675 Brownian bridge simulations within the 
accept-reject step (iii) at each iteration of the algorithm. By comparison RWPF1 and 
RWPF2 simulated respectively 910 and 1000 particles and required on average 1025 and 
850 Brownian bridge simulations to generate the random weights at each iteration. 

Note that the comparative CPU cost of the four algorithms, and in particular that 
of EPPF and ESPF as compared to RWPF1-2 depends on the underlying diffusion path. 

k ■ ■ 

The acceptance probabilities within EPPF and ESPF depend on the values of x^' 3 and 
Xt i+1 , and get small when both these values are close to 0(mod 2tt). (in the long run the 
diffusion will visit these regions infrequently and will stay there for short periods.) Thus, 
simulated paths which spent more (or less) time in this region of the state-space would 
result in EPPF and ESPF having a larger (respectively smaller) CPU cost. 

We compared the four filters based on the variability of estimates of the mean of the 
filtering distribution of the state across 500 independent runs of each filter. Results are 
given in Figure El while output from one run of RWPF2 is shown in Figure [U The com- 
parative results in Figure El are for estimating the mean of the filtering distribution at each 
iteration (similar results were obtained for various quantiles of the filtering distribution). 
They show RWPF2 performing best with an average efficiency gain of 15% over RWPF1, 
50% over ESPF and 200% over EPPF. Interpretation of these results suggest that (for 
example) ESPF would be required to run with N = 750 (taking 1.5 times the CPU cost 
for this data set) to obtain comparable accuracy with RWPF2. 

Varying the parameters of the model and implementation of the algorithms will affect 
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the relative performance of the algorithms. In particular increasing (or decreasing) a 2 , 
the variance of the measurement error, will increase (respectively decrease) the relative 
efficiency of EPPF relative to the other filters. Similar results occur as A, is decreased 
(respectively increased). The relative performance of the other three algorithms appears to 

(k) (k) 

be more robust to such changes. We considered implementing EPPF with fll = w\ ■ ; and 
also using an Euler rather than an Ozaki approximation of the sine diffusion to construct 
the proposal distribution for RWPF1-2, but neither of these changes had any noticeable 
effect on the performance of the methods. We also considered re-sampling less often, setting 
C = N/4 in step PF1 of the RWPF algorithm (so re-sampling when the effective sample 
size of the $ j) s was less than N/4) and this reduced the performance of the algorithms 
substantially (by a factor of 2 for RWPF1-2). 

We also investigated the effect of increasing the amount of time, A, between observa- 
tions. To do this we used the above data taking (i) every 10th; or (ii) every 20th time-point. 

To measure the perfor mance of the filter for t hese different scenarios we used the Effec- 
tive Sample Size (ESS) of Carpenter et al. ( 19991 ). ESS is calculated based on the variance 
of estimates of posterior means across independent runs of the filter, but this variance is 
compared to the posterior variance to give some measure of how many independent draws 
from the posterior would produce estimators of the same level of accuracy. We focus on 
estimates of the posterior mean of the state at observation times; and if s 2 is the sample 
variance of the particle filter's estimate of E[X t . across 100 independent runs, and a 2 
is an estimate of Var[X ti \yi : i], then the ESS is a 2 /s 2 . Note that comparing filters by their 
ESS is equivalent to comparing filters based on the variance of the estimators. 

Table [3] gives ESS values for the different values of A. We see that the ESS values 
drops dramatically as A increases, and the filter is inefficient for A = 20. This drop 
in performance is due to the large variability of the random weights in this case. The 
variability of these weights is due to (a) the variability of 

exp(- t +1 g(W s )ds\, (26) 



across different diffusion paths; and (b) the Monte Carlo variability in estimating this for a 
given path. To evaluate what amount is due to (a), we tried a particle filter that estimates 
(T26]) numerically by simulating the Brownian Bridge at a set of discrete time points (for 
this example we sampled values every 1/2 time unit) and then using these to numerically 
ev aluate the integral. This app r oach is closely r elated to the importance sampling approach 
of Durham and Gal lant (120021 ): iNicolaul <j2002h . see Section O The results for this filter 
are also given in Table [3] (note the ESS values ignore any bias introduced through this 
numerical approximation), and we again see small ESS values, particularly for A = 20. 
This filter's performance is very similar to the RWPF, which suggests that the Monte Carlo 
variability in (b) is a small contributor to the poor performance of the RWPF in this case. 

Finally we tried introducing pseudo observations at all integer time-intervals where cur- 
rently no observation is made. The RWPF is then run as above, but with no likelihood 
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Filter 


A = 10 


A = 20 


RWPF2 


73 


5 


Discretisation 


80 


12 


pseudoRWPF2 


923 


933 



Table 3: Comparison of filter's mean ESS values for different time intervals between obser- 
vations (A). Results are for the Random Weight Particle Filter using GPE-2 (RWPF2), a 
filter that numerically approximates the weight through discretising the diffusion process 
(Discretisation), and the RWPF after introducing uninformative observations at unit time 
intervals (pseudoRWPF2). 



contribution to the weight at the time-points where there are these uninformative obser- 
vations. The idea is that now A = 1, so that the variance of the random weights is well- 
behaved, but we still have adaptation of the path of the diffusion in the unit time-interval 
prior to an observation to take account of the information in that observation. Results are 
again shown in Table [3], and the ESS values are very high (and close to the optimal value, 
that of the number of particles, 1000). Note that the computational cost is only roughly 
doubled by adding these extra pseudo observations; as the total computational cost for 
the simulation of the Brownian bridge is unchanged. These results are reasonably robust 
to the choice of how frequently to introduce these uninformative observations (results not 
shown) . 

5.2 Analysis of the Cox process 

We now consider applying the random weight particle filter (RWPF) to Example 2 from 
Section [2], the OU-driven Cox process. The data we analysed is given in Figure [3](top). 
It is either impossible or difficult to adapt the other two EA-based particle filters (the 
EPPF and the ASPF) to this problem. For instance we cannot implement EPPF as the 
likelihood function is not tractable. As such we just focus on the efficiency of the RWPF 
in estimating the filtering distribution of \X t \. 

Our implementation of the RWPF was based on proposing particles from the prior 
distribution, so = w\ and q{x ti+1 |xf , Vi+x) is just the OU transition density p{x ti+l \x{). 
We simulated the random weights by GPE-2. We calculated the filtering density at each 
observation time, and also at 56 pseudo-observation times chosen so that the maximum 
time difference between two consecutive times for which we calculated the filtering density 
was 0.1. This was necessary to avoid the number of Brownian bridge simulations required 
to simulate the weights being too large for long inter-observation times, and also to control 
the variance of the random weights (see above). The likelihood function for these non- 
observation times is obtained by removing v{x ti ) from (jl]). 

We set the number of particles to 1, 000 and resampled when the ESS of the f3^s was 
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Figure 3: Top: Simulation from the Cox process of Example 2 and results from anal- 
ysis by the RWPF. The path of the absolute of the underlying diffusion (black line); 
observed arrival times (green dashes); filter estimates from the RWPF (red circles); and 
90% credible interval for the absolute of the diffusion (red dashed line). (Whilst for clarity 
they are shown for all time, the credible intervals were only calculated at and apply for 
times where filter estimates are shown.) Bottom: ESS of the RWPFs weights (defined as 

(J2f=i w i^) 2 / J2j=i( w i^) 2 ) over time. The dramatic increases in the effective sample sizes 
correspond to re-sampling times. 



less then 100 (C = iV/10 in step PF1 of the algorithm in Section [3]). Whilst results for the 
sine diffusion suggest that this will result in an algorithm that re-samples too infrequently, 
we chose to have a low threshold so that we could monitor the performance of the particle 
filter by how the ESS of the particle filter weights decay over time. The results of one run 
of this filter are shown in Figure Eltop). The computational efficiency of this method can 
be gauged by Figure [3] (bottom) where the ESS of the w^'s is plotted over time. 

6 Discussion 

We have described how recent methods for the exact simulation of diffusions and the 
unbiased estimation of diffusion exponential functionals can be used within particle filters, 
so that the resulting particle filters avoid the need for time-discretisation. Among the 
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approaches we have introduced special attention was given to RWPF which implements 
an auxiliary particle filter, but simulates the weights that are allocated to each particle. 
We showed that this methodology is equivalent to an auxiliary particle filter applied to 
appropriately expanded model. We expect that this methodology will have interesting 
applications to different models than those considered in this paper, which however involve 
intractable dynamics or likelihoods. 

We have focused on the filtering problem, estimating the current state given observa- 
tions to date. However, extensions to prediction are trivial - merely requiring the ability 
to simulate from the state equati on, which is pos sible via the EA algorithms. It is also 
straightforward to use the idea of iKitagawal (119961 ). where each particle stores the history 
of its trajectory, to get approximations of the smoothing density (the density of the state 
at some time in the past given the observations to date). 

Note that while particles store values of the state only for each observation time, it 
is straightforward to fill in the diffusion paths between these times to produce inferences 
about the state at any time. A particle approximation to the distribution of (X s ,tj_i < 
s < ti), conditionally on the data can be constructed using the current set of weighted 



particles {(x-i\, itf )}f = i with weights {w^ 1 }, as follows. Firstly we need to introduce 

some notation; we denote by the value of the particle at time t^_i from which the jth 
particle at time ti is descended. The particle approximation is given by a set of weighted 



,(?) 



paths {(x s ,ij_i < s < tj)^}^! with weights {w^}. Each path is a diffusion bridge 



starting from x ^-,^. and finis hing at x) J> a nd it can be simulated using EA, as described 
Beskos et al.l (j2006al ) and iBeskos et al.l (j2005bl ). In observation regimes (A) and (B) 



in 



the EA is applied to simulate a diffusion bridge with density w.r.t. the Brownian bridge 
measure given by exp{— f*' ^ 0(X s )ds}, whereas in regime (C) the corresponding density 

is exp{— f * (0(X S ) + z/(X s ))ds}. This representation can be directly exploited to draw 
inferences for any function of a finite skeleton of X in-between observation times. 



Appendix A: The layered Brownian motion 



The algorithm proposed in lBeskos et al.l (l2005bl ) starts by creating a partition of the sample 
space of W for the given W = x and W t = y. Writing x = (x\, . . . , x^), for a user-specified 
constant a > a sequence of subsets of R d is formed as Aj = {u = (u±, . . . , «d) ■ 

min(x ti , yi) — ja < Ui < max(x ti , y,i) + ja}, j > 0, where UjAj = R d . This sequence defines 
a partition of the sample space of the form UjZiDj, where a path belongs to Dj if and only 
if the path has exceeded the bo unds determined by A,_x but not the bounds determined 
by Aj. In IBeskos et al.l (l2005bl ) it is shown how to simulate the random variable which 
determines which of the Djs W belongs to, and how to simulate W at any collection of 
times conditional on this random variable, the layered Brownian bridge construction. Since 
g is assumed continuous, knowing W G Dj is sufficient to determine Uw and Lw which 
satisfy ( |T6l) . In fact, in the simplified setting where g is bounded, as in the sine diffusion 
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of Example 1, the layered Brownian bridge construction can be avoided since it is easy to 
choose f/w and L w independently of W. 



Appendix B: Proof of Theorem Q] 



t K 



k\p(k I U W ,L W ] 



Then, ( TH1 is established as follows: 
E[/ 2 | U w , L w ] = E[ E[/ 2 | k, W] ] = E 
= E 

oc 

= £ 



i2fi 



[k\p(k I U W ,L W )) 2 \J 



(u w - g (w s )y 



ds 



t K 



[K\p(k | U W ,L W )} 



(U w -g(W s )) 2 ds) \U W , 



k=0 



p(k | U-w, L w )k\ 2 



E 



j\u w -g(W s )) 2 ds^J \U W ,L W 



Fubini's theorem and dominated convergence are used above (valid since the integrands 
are positive a.s.). ( !22l) is obtained using the following result (which can be easily proved 
using Jensen's inequality). Let /j > for i = 1, 2, . . .. Then the sequence of piS which 
minimize Yli^o fi/Pi under the constraint J^pj = 1 is given by pi = \fYil Yl VJi- 

Appendix C: Proof of Theorem [2] 

GPE-1< e _iw< so that the result holds if E[e _Lw '] < oo, where the expectation is w.r.t. a 
(^-dimensional Brownian bridge from x at time to y at time t. However 



POO 

- / P[e~ Lw > w}dw 
Jo 

poo poo d 

= / P[L W < -\ogw]dw < / P[5 V(l + Mi) > \ogw]dw 
Jo Jo i=1 

where = sup 0<s<t \ Wi\ using the the growth bound in Theorem [2 Furthermore, 

poo d poo d 

/ P[6 V(l + Mi) > \ogw]dw < / VP[5(l + Mi) > d~ 1 logw]dw 
Jo ' j=1 Jo i=1 

poo d, 

= / VP[Af, > (dd)- 1 \ogw - l]dw . 
Jo i=i 
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It remains therefore to bound the d integrals on the right hand side of this expression. How- 
ever from the Bachelier-Levy formula for hitting times for Brownian motion and bridges, 

P[Mj > v] < exp {-2(v - max{xi, Vi}) 2 /t} + exp {-2(min{xj, y { + v}) 2 /t] 

and so 

P[Mi > (d5)- 1 logw - 1] < exp{-2((d5)- 1 logu; - 1) - max{x h yi}) 2 /t} 

+ exp{-2(mm{x hyi + (dS)' 1 log w - l)}f/t} 

which recedes like w~ klogw as w — > oo thus concluding the proof. 

Appendix D: EPPF and ESPF for Example 1 

EPPF generates the new particles according to the following procedure: 

(i) choose one of the current particles x\ k , where particle j is chosen w.p. j3^\ 

(ii) propose x tl+1 from Normal with mean x < f , '' 3 \ and variance A,; 

(hi) accept this proposal w.p. exp(— cos(x^ +1 ) — 1); if proposal is rejected return to (i). 

(iv) accept this proposal with probability 



E 



A, 

exp { - / <P(W s )ds 



o 



(k- ■) 

where expectation is with respect to the law of a Brownian Bridge from Wq = x\ 
and = Xt l+1 , and (ft is given in (J7J). If the proposal is rejected return to (i), 
otherwise x ti+1 is the new particle at time t i+ i with weight w i+1 = f(yi + i\x i+1 ). 



(iv) is performed using retrospective sampling as described in Beskos et al. ( 2006al ). 



ESPF proceeds as above but with steps (i) and (ii) replaced by the step 
(i') propose (x[ kl ' 3 \ x tl+1 ) according to the density proportional to 

(/,,,., yV^- x t l ' j) ) [ f (.o , - //) ;: 



exp < — cos{x i ) ' — > exp 



2(cr 2 + (Aj)) \ *\ 2r 

where 7] = (a 2 + A^' 1 (x^ a 2 + A^ +1 ), and r = a 3 At/ (a 3 + A,). 
The algorithm is repeated until iV values for x ti+1 are accepted, each with weight 1/JV. 
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Appendix E: Proposal Distribution for Example 1 

Consider a diffusion satisfying SDE (CQ), with d = 1 for simplicity. The Ozaki approximation 
of this SDE is based on a first order Taylor expansion of the drift about some value x. For 
the sine diffusion of Example 1, we get the following approximating SDE 

dX s = — cos(x)[x — tan(a;) — X s ] ds + dB s . 

So X s — (x — tan(x)) is an OU process as defined in Example 2 with p = cos(x) and a = 1. 
To calculate q{x ti+1 \xf \ yi+i) we compute the product of the transition density given by 
the Ozaki approximation about x = xf' and the likelihood function f (y i+ i\x ti+1 ) . Defining 
r 2 = (l-exp{-2cos(x 4 (i) )A i })/(2cos(a; l (i) )), andr] = a;S J ' ) -tan(x?' ) )(l-exp{- cos(4 j) )Ai}) 
we get that q{x ti+1 \x^\ Ui+i) is Normal with mean (77a 2 + yi + \T 2 ) /(r 2 a 2 ) and variance 
rfr 2 / (if + t 2 ). Furthermore we calculate oc w^J\f T 2 +a 2 (yi+i — rf). 



Appendix F: Central Limit Theorem 

For notational simp licity, we consid er a special case of our particle filter, chosen to resemble 



those considered in IChopinl (12004] ) . We choose our proposal density for time U + \ to have 



0i = y ij'i an d we assume iid sampling of A" t in step PF1. The particle filter of iChopin 



(120041 ) splits up simulating particles at time t i+ % into (i) a resampling of particles at time t 
and (ii) a propagation of each of these particles to ti me . Our a ssumption of iid sampling 
is equivalent to the multinomial resampling case of Chopin ( 2004I ). ( The conditions for th e 



central limit theorem are the same if the residual sampling methods of lLiu and Chenl (Il998h 
but the variances differ.) For simplicity we consider observation model (A) or (B), though 
the result extends easily to observation model (C). 

Let ey> = ( 

x ti\ x u-i^)> "where fcjj is the index sampled in step PF1 when simulating 
the j particle at time U and 9^ is the jth particle at time U together with the particle 
at time from which it is descended. Also let denote conditional expectation given 
6>j. Similarly, let Hi{6i) = x i: tj), and denote by Ri the unbiased estimator of 

Pi(9i), i.e. E[i2j] = Hi{9i). An important quantity is &f(9i) = Var(_Rj). 

We define Ejk/>] and Vaii((p) to be the posterior mean and variance of an arbitrary 
function ip{6) at time i, and consider Particle Filter estimates of Ej [<p]. Let fci(9i) be the 
density p(x ii _ 1 \y\.i-i)q{x ti \x ti _ 1 ). Finally define E 9i [<p] and Var 9 .(yj) to be shorthand for 
the conditional expectation and variance of <p(9i) with respect to q(x t .\x ti _ 1 ) (which are 
functions of x 4i l ). We denote || ■ || to be the Euclidean norm and define recursively $; to 
be the set of measurable functions ip such that for some 5 > E # .[||/i;ii:^|| 2+<5 ] < 00, and 
that the function Xt i _ 1 i— > E 9i [hi[iitp\ is in 

Theorem 3. Consider a function ip; define Vq = Var^^^Lp) , and by induction: 

Vi{tp) = Vt-i {EM} + Ei-i { Var qi {cp)} , for 1 > 0, (27) 
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Viifp) 



for i > 0, 



V i (i f ) = V i (<P)+Var l (<f ) ), fori>0. 

1 belongs to (C2) E^[hja 2 ] < oo; and (C3) E^ [a i ^h i ] 2+s 



(28) 
(29) 



Then if for alii (CI) x u h-> 1 belongs to <Pj/ (G'^y It^J^cxfj < oo; and (U3) 11%. [o^AijJ < 
oo /or some 5 > 0; inen /or any psfj, Ej[<£>] and are /imte and we /iai>e iae following 
convergence in distribution as the number of particles, N , tends to infinity: 



N l/2 



N 



W 



N 



W 



U) 



Comment Equations (|27|) - (|29p refer to the changes in variance of the weighted particles 
due to the propagation, weig hting an d resam pling stages at iteration i. Only ( !28l) differs 
from the respective result in IChopinl (12004] ) . and this is due to the second term on the 
right-hand side, which represents t he increase in variance due to the randomness of the 
weights. Condition CI is taken from lChopinl (120041 ) and applies to standard particle filters; 
conditions C2 and C3 are new and are conditions bounding the variance of the random 
weights which ensures that Vi(p) is finite. 



Proof. We adapt the induction proof in lChopinl (120041 ) . considering in turn the propagation, 
weighting and resampling steps Our filter differs from the standard particle filter only in 
te rms of the weig hting step; and therefore we need only to adapt the result of Lemma A2 
in IChopinl (120041) . In fact, (127)1 and (I29*|) are identical to the corresponding quantities in 
Chopi nl (120041 ). therefore it remains to show (I2"8~l) . We define the constant K = E^Rihi] 



and ip* = Rjhj( (fi — Ej((p))/K. Within the enlarged signal space framework, we can apply 
Equation (4) of lChopinl (I2004h . to give: 



V 4 (<p) = ^(^) = ^-i{E g J^]} + E l - 1 {Var, i (^)}. 

Now we can calculate E qi [95*] by first taking expection over the auxiliary variables (condi- 
tional on 6j). This gives E ? .[y>*] = E gi [/ij/tj(</? — Ei(ip)) / K]. Similarly we get 

Var 9i (^*) = Var gi (E[^(^-^(^))/K])+E ft (Var[i?A(^-^M)/^]) (30) 
= Var,,^^ - Ei((p))/K) + E qi (a 2 h 2 (cp - E^f/K 2 ). (31) 

(Here the expectation and variance in ( 1301) are w.r.t. the auxiliary variables). Combining 
these results gives ( 1281) . The regularity conditions (CI) - (C3) translate directly also. □ 
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